Jammed Spheres: Minkowski Tensors Reveal Onset of Local Crystallinity 
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The local structure of disordered jammed packings of monodisperse spheres without friction, 
generated by the Lubachevsky-Stillinger algorithm, is studied for packing fractions above and below 
64%. The structural similarity of the particle environments to fee or hep crystalline packings (local 
crystallinity) is quantified by order metrics based on rank-four Minkowski tensors. We find a critical 
packing fraction <j> c « 0.649, distinctly higher than previously reported values for the contested 
random close packing limit. At <j> c , the probability of finding local crystalline configurations first 
becomes finite and, for larger packing fractions, increases by several orders of magnitude. This 
provides quantitative evidence of an abrupt onset of local crystallinity at <p c . We demonstrate that 
the identification of local crystallinity by the frequently used local bond-orientational order metric qe 
produces false positives, and thus conceals the abrupt onset of local crystallinity. Since the critical 
packing fraction is significantly above results from mean-field analysis of the mechanical contacts for 
frictionless spheres, it is suggested that dynamic arrest due to isostaticity and the alleged geometric 
phase transition in the Edwards framework may be disconnected phenomena. 
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PACS numbers: 05.20.Jj statistical mechanics; 61.20.-p structure of liquids; 45.70.-n granular systems 



Now classic experiments showed that disordered sphere 
packings can only be prepared up to a maximum packing 
fraction of </>rcp ~ 0.64 [H, Q. This packing fraction, re- 
ferred to as random close packing (RCP), is substantially 
lower than the packing fraction <fif cc w 0.74 of the dens- 
est crystalline sphere packing. Maximal packing fractions 
close to </>rcp have been shown for several experimental 
protocols [3[; protocols inducing local crystallization are 
able to reach higher packing fractions [1] . Numerical pro- 
tocols to generate static sphere packings both below and 
above 0.64 are the Lubachevsky-Stillinger (LS) algorithm 
@ and the Jodrey-Tory algorithm @. 

The nature and the existence of a transition near </>rcp 
are disputed. As sphere configurations with packing frac- 
tions between </>rcp an d </>f C c evidently exist |7j , an alleged 
structural transition must be due to either a vanishing 
configuration space density of these states or to the in- 
ability to reach these within the considered ensemble or 
by the given dynamics. Within the framework of equi- 
librium (thermal) hard spheres, the concept of RCP has 
been related to the terminus of the branch of metastablc 
states avoiding crystallization; divergence of pressure is 
reported to occur at <f> = 0.640±0.006 g| and « 0.65 @. 
By contrast, in an athermal statistical ensemble where 
the role of energy is played by volume [l(J, a mean- field 
study based on mechanical contact numbers has reported 
0r,CP ~ 0.634 An order/disorder transition in an 



athermal ensemble has been demonstrated for a lattice 
model [l2| . Support for the phase transition scenario is 
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deduced from the fact that the volume fraction of polytc- 
trahedra increases with packing fraction up to (f> sw 0.646 
and then decreases, as these structures transform into 
crystalline order [l3|, EI • 

In addition to these phase transition scenarios, where 
0rcp is interpreted as the density of the disordered 
phase at coexistence, the notion of the maximally ran- 
dom jammed state (MRJ, [lj|) has been proposed (as the 
maximally disordered state amongst all jammed pack- 
ings); with respect to a number of common measures of 
order, the MRJ packing fraction has been estimated as 
0MRJ « 0.63 0. 

The resolution of the RCP problem relies on suitably 
defined order metrics to quantify packing structure. A 
common approach to local structure characterization is 
by analysis of nearest neighborhoods [III [I?} ■ Often, the 
bond-orientational order metrics qi defined by Steinhardt 
et al. [HI are applied to sphere packings However, 
these and other neighborhood-based order metrics (l9| 
have shortcomings. First, they suffer from the ambiguous 
definition of the nearest neighborhood [2(J. Second, in 
their common use as single scalar order metrics 0, [H| , 
they are insufficient to conclusively distinguish order and 
disorder [2l[ . A non- negligible fraction of non-crystalline 
environments is often incorrectly identified as crystalline 
(false positives), since their q§ are close or identical to 
the reference qe values in crystals. 

Alternatively, the structure of monodisperse sphere 
packings can be characterized by analysis of the Voronoi 
cells of the particle centers, see Fig. [TJ Suitable morpho- 
logical descriptors, such as Minkowski tensors [22l. |23|. 
can then be used to quantify the cell shape and hence 
the local structure. Here, we show that crystalline or- 
der metrics can be constructed from rank-four Minkowski 
tensors of the Voronoi cells that give stringent criteria for 
fee or hep crystalline order. For jammed sphere packings 
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FIG. 1. (Color online) (left) The Voronoi diagram of a 
packing associates convex cells with the individual particles. 
Each Voronoi cell K contains one particle k such that the 
distance of any point p £ K to particle k is smaller than the 
distance of p to any other particle, (right) For the evaluation 
of qe, we define the set of nearest neighbors of particle k as 
those 12 other particles (6 in 2D) which are closest to k. 



generated by the LS algorithm, these order metrics re- 
veal an abrupt onset of crystallinity at a critical packing 
fraction <fi c w 0.649. 

Eigenvalue ratios of rank-two Minkowski tensors quan- 
tify anisotropy of the particle environments in jammed 
bead packs [22| . The Voronoi cells of the fee and hep 
close packing arc isotropic, in the terminology of Ref. [22J, 
while cells found in disordered packings typically are not. 
However, rank-two tensors are insufficient to distinguish 
different types of isotropic cells. These cell types can 
be classified using the rank-four Minkowski tensor W^ 0,4 . 
The tensor W®' of a Voronoi cell is given as the sum 
of tensor products of the facet normals, weighted by the 
facet areas, 
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where rii :~ (n(/))i with i — 1,2,3 are the cartesian 
components of the facet normal, and a(f) is the surface 
area of facet /; further, A := a(f) is the total sur- 
face area; all sums run over the facets of the Voronoi cell 
K. In close analogy to the stiffness tensor of continuum 
mechanics, symmetry under permutations of indices al- 
lows the reduction of W®' A to a 6 x 6 symmetric matrix 
[24| . The six eigenvalues (ft, . . . , of this matrix are 
dimcnsionless due to normalization by A -1 and are rota- 
tional invariants pBT ] . A concise quantitative measure for 
the similarity of a given Voronoi cell K to the Voronoi 
cell Kf cc of a crystalline fee packing is given by the fee 
crystalline order metric 
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An analogous order metric Ah cp is denned for hep cells. 

Figure [U supports our claim that Af cc and Ah cp mea- 
sure deviations of a Voronoi cell's shape from the ideal 
fee or hep cell. The vertices of an ideal fee or hep lattice 
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FIG. 2. (Color online) The relationship between root mean 
square particle displacement from an ideal fee (red solid line) 
and hep (green dashed) lattice and and the corresponding val- 
ues of Af cc and Ahcp respectively is approximately linear. Er- 
ror bars represent standard deviations among cells due to the 
statistical displacement of the particles. The dashed curves 
above and below the data are linear functions with slope 1 
and 0.1 respectively. 



are displaced by small random vectors. Figure [2] shows 
averages and standard deviations (as error bars) of Af cc 
and Ahcp as function of the root mean square displace- 
ment (RMSD) from the ideal lattice points, demonstrat- 
ing an approximately linear relationship between devia- 
tions from the ideal crystalline shapes and the crystalline 
order metrics Af cc and Ahcp- The quantitative agreement 
between the functions for hep and fee cells justifies the 
use of the same threshold value for selecting both fec- 
and hep-like cells from a packing. 

The crystalline order metrics Af cc , Ah cp are used to 
identify crystalline clusters in jammed sphere packings 
generated by the Lubachevsky-Stillinger protocol [ff. 
Figure [3] shows, as key result of this paper, that (a) crys- 
talline fee and hep order is absent for packing fractions 
below a critical value <f> c « 0.649; and that (b) above 
(j) c the fraction of crystalline fee or hep in LS simula- 
tions is non-zero and rapidly increases by several orders 
of magnitude. As expected for an athermal system, no 
systematic difference between the number of hep and fee 
cells is observed, in contrast to crystallization dynamics 
in equilibrium [26J. 

We measure the fraction of fee cells as n{ cc (cf>) = 
[Af cc ]o.5/A, where iVf cc is the number of cells with Af cc < 
8 = 0.005, in a sample of TV = 4 x 10 4 spheres. [Af cc ]o.5 
is the median of Af cc over M ss 20 packings of similar 
(j), see Fig. [3] (In general, for a random variable X with 
a probability density /(A), the symbol [X] p denotes the 
p-quantile, i.e. the value F~ 1 {p) with the cumulative dis- 
tribution function F{X) = J* /(£)d£.) The accuracy 
of our results is, for small rif cc and n\ lcp , limited by the 
finite system size, preventing the measurement of prob- 
abilities smaller than 1/A. (We further note, that the 
number of isotropic cells, in the terminology of Ref. [22j], 
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FIG. 3. (Color online) Fractions nf cc (<^) (red) and nh cp (0) 
(green) of fee- and hep-like cells as function of the packing 
fraction <f> (cell types identified by W-f' 4 ). Below 4> c ~ 0.649, 
crystallinity is negligeable, and was found, amongst 3000 sim- 
ulations with cj) G [0.56,0c]) with rates substantially smaller 
than the inverse system size. Above <f> c , the probability 
of crystalline cells increases by orders of magnitude. Each 
data point is computed from M ~ 20 packings of 4 x 10 4 
spheres each. Horizontal error bars correspond to the vari- 
ations in <f) encountered for the same growth rate 7 of the 
LS algorithm vertical error bars represent the interval 

([JVfc C ]0.25,[^Vfcc]0.75). 



FIG. 4. (Color online) The fee crystalline order metrics Af cc 
of the most fcc-like cells in each packing, quantified by the 
first percentile [Af cc ]a; with x — 0.01 of the Af cc distribu- 
tion in each individual packing. (Each datapoint represents 
one packing). The insert shows variations of the estimate for 
4> c obtained as the intersection point of fitted straight lines, 
when x is varied. The black bars indicate published estimates 
of (/>rcp: (a) Anikeenko and Medvedev's analysis of tetrahe- 
dral configurations [l3|, [l4| , (b) Bernal's analysis of steel ball 
bearings Q], and (c) the contact number analysis by Song et 

al. m. 



vanishes below <p c , and becomes nonzero at <j) c .) 

The values of rif cc and Uh cp depend, of course, on the 
choice of the threshold 8. The increase in local crys- 
tallinity is, however, also evident in the lowest occur- 
ring values of Af cc . Figure |4] shows the first percentile 
[Af cc ]o.oi as a robust estimate for the most crystal-like 
cell, i.e. the lowest occuring value of Af cc . A sharp drop 
of [Af cc ] 0.01 is observed for <j) > <fi c ; in packings below 
</> c , the most fcc-like cells are substantially different from 
fee cells, while above 4> c the differences quickly decay 
to close to zero. The value of 4> c (x) is estimated by 
the intersection of two straight lines fitted to the data 
for [Afcc] x- The insert of Fig. [4] shows the (j> c estimates 
extracted by this approach, giving C G [0.6492, 0.6499] 
for x E [0.001,0.1]. These values of <f) c are larger than 
published values for the RCP or the MRJ packing frac- 
tion. Importantly, the data of Fig. [4] demonstrate that 
the drastic increase in nf cc in Fig. [3] is a robust result 
that is not sensitive to the value of the threshold 8. The 
value of 8 is, within bounds, an irrelevant parameter. We 
do not observe differences between packings of N = 10 4 
and N = 4 x 10 4 particles besides improved statistics. 

The observed abrupt appearance of crystalline cells at 
4> c is difficult to observe using the bond-orientational or- 
der metrics qi and wi developed in seminal work by Stein- 
hardt et al. [161 ] . Most frequently, q e is considered, which 
is deemed particularly sensitive to formation of fee; it is 
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with the spherical harmonics Y/ m , the polar angles 9^ 
and tfjk of the bond vector between particles j and k, 
and (•) denoting the average over the 12 closest neighbors 
j of k (Fig. EE]). 

Figure [5] shows probability distributions of q§ values 
observed in LS packings above 4> c that demonstrate the 
principal deficiency of using only q$ as an order metric. 
The frequency /((ft) of qe values develops sharp peaks at 
the values corresponding to fee (q§ — 0.57452) and hep 
(<76 = 0.48476) for <fi > <p c , not present for samples with 
(f> < (f) c . These peaks, however, sit on top of a dominant 
background of non-crystalline cells. The data clearly 
show that cells (false positives) exist which are distinctly 
different from fee but that are identified as fee by q^, 
i.e. \qe — q f 6 cc \ < 5 • 10 -4 . For example, the cell displayed 
in (d) has eleven facets, several of which are five-sided; 
analogous hep examples exist. If cells that are identified 
as either fee or hep by W®'^ are excluded from the qe dis- 
tribution, these peaks vanish; the residual smooth distri- 
bution represents the non-crystalline background. Thus, 
for reliable detection of crystallinity, more information is 
required than contained in qe alone. This can be achieved 
by using multiple qi metrics and their distributional prop- 
erties [21], or more specialized bond-orientational order 
metrics such as 6* fcc , 8 hcp [T3|. Minkowski tensors, such 
as used in the present study, represent a more general 
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FIG. 5. (Color online) Abundances of all cells with a specific 
Q6 value (/ail, red solid curve) and of the subset of cells that 
are very clearly neither fee nor hep, i.e. with Af cc > 0.015 
and Ahcp > 0.015 (/f p , blue dotted curve). The values of the 
blue, dotted curve at q^ , Qg Cp are finite even though genuine 
fee and hep cells are excluded. This clearly demonstrates 
that qa produces false positives (fp), that is cells that are not 
crystalline but identified as such by q$. Only the difference 
fa := /ail — /f P consists of truly crystalline cells. The cells 
depicted represent (a) an ideal hep cell, (b) an ideal fee cell, 
and cells identified by q@ , but not by Af cc or A hcp , as (c) hep 
and (d) fee. The data is averaged over ten configurations, each 
consisting of N = 4 x 10 4 spherical particles, with packing 
fractions 4> m the interval [0.656,0.660], well above C . See 
Fig. [T] for the definition of particle neighborhood. 



approach; it is not necessary to choose a set of neighbors 
or bonds associated with a particle in order to evalu- 
ate the Minkowski tensors, and the Minkowski tensors 
are continuous functions of the particle coordinates. At 
the same time, they contain the information necessary 
to discriminate in a robust and specific way between dis- 
ordered structure and different types of crystallinity. It 
can be shown that a relation exists linking the rotational 
invariants of higher-rank Minkowski tensors to variants 
of the bond-orientational order metrics qi , wi which have 
been amended by weighting factors proportional to the 
Voronoi facet areas [20| . 

In conclusion, we have demonstrated that local crys- 
tallinity in Lubachevsky-Stillinger sphere packings sets in 
when the packing is compactified beyond a critical pack- 
ing fraction </> c w 0.649. The packing fraction </> c marks 
the density below which LS configurations show no de- 



tectable degree of local crystallinity. Compactified above 
this limit, the system responds by the formation of local 
crystallinity. 

The value <fi c ~ 0.649 is higher than experimental es- 
timates for the RCP limit [j], @] and than the predic- 
tion based on mechanical contact numbers (Tlj . but also 
than the MRJ packing fraction [p|- However, <p c is close 
to the packing fraction of 0.646 where polytetrahe- 
dral aggregates are most prevalent (Fig. 3 in [13|); the 
crystalline order metrics thus identify the conversion of 
polytctrahedral aggregates into crystalline structure, de- 
tected indirectly by Refs. [13, 13- The small but signif- 
icant discrepancy between the critical packing fractions 
< 0.64 of Ref. [J 0, H [ll| on the one hand and « 0.65 of 
Refs. [13, EH and of the present work on the other raise 
the caution that the mechanisms of dynamic arrest and 
isostaticity may be distinct from the alleged geometric 
order-disorder transition. 

Given the nonequilibrium nature of jammed matter, 
one might be tempted to attribute the observed differ- 
ence in packing fraction to details of the preparation 
protocol. The critical packing fraction of « 0.65 is, how- 
ever, not specific to the LS algorithm. For example, for 
packings generated using the force-balancing 'split algo- 
rithm', the geometric (rather than mechanical) contact 
number exhibits an anomaly close to <j) ~ 0.65 (Fig. 14 
of Ref. [27j). Data by Bargiel and Tory for Jodrey-Tory 
packings can be successfully fitted with a critical packing 
fraction of « 0.6495 [13] ■ Recent results by Klumov et al. 
for Jodrey-Tory and Lubachevsky-Stillinger packings, in 
terms of quantiles of the w§ distribution, fix the geomet- 
ric transition around 0.65 [2l[, but do not exclude 
crystallinity below <j> c . 

Future work needs to focus on the precise nature of 
the geometric transition occurring at <p c . Is the first- 
order phase transition scenario viable, and if so, what are 
the coexistence densities? What is the signature of the 
transition in the Af cc distribution? How does the local 
structure (quantified by Minkowski tensors) relate to the 
observed "Kauzmann" density (Fig. 8 of Ref. flij)? 
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